Also known at the Ordinary Least Squares (OLS) model
It aims to fit a regression line which minimises the squared differences between observed and predicted values
Equation of the line
The equation of a linear model (lm) is given by:
\[
y_i = \beta_0 + \beta_1 x_i + \epsilon_i
\]
where:
\(y_i\) is the predicted value of the response variable.
\(\beta_0\) is the intercept.
\(\beta_1\) is the slope coefficient, representing the change in the response variable for a one-unit change in the explanatory variable.
\(x_i\) is the value of the explanatory variable.
\(\epsilon_i\) represents the residuals of the model
Limitations of general linear models
Q. What are the assumptions of a general linear model?
Assumes a linear relationship
Assumes normal distribution of the residuals
Assumes homogeneity of the residuals
Independence - assumes observations are independent of each other
Normal distribution of the residuals
The linear regression line is the most likely line given your data if we assume each data point comes from a hypothetical bell curve centered on the regression line
Normal distribution
Normal distribution of the residuals
When using our residuals to calculate standard error, confidence intervals and statistical significance these are assumed to be drawn from a normal distribution with mean zero and a constant variance.
This implies that the residuals, or the distances between the observed values and the values predicted by the linear model, can be modeled by drawing random values from a normal distribution.
Which literally means that \(y_i\) is drawn from a normal distribution with parameters \(\mu\) (which depends on \(x_i\)) and \(\sigma\)(which has the same value for all measures of \(Y\)).
Real data
How often do these assumptions really hold true?
A good fit:
Real data
How often do these assumptions really hold true?
An ok fit:
Real data
How often do these assumptions really hold true?
A poor fit:
Assumption testing
Residuals vs Fitted
Purpose: To check for linearity and homoscedasticity.
Interpretation: Look for a random scatter of points around the horizontal line at y = 0. A funnel-shaped pattern suggests heteroscedasticity.
QQ-plot
Purpose: To assess the normality of the residuals
Interpretation: Points should fall along the diagonal line. Deviation from the line indicates non-normality of residuals.
Scale-location
Purpose: To check homoscedasticity and identify outliers
OK: Error variance appears to be homoscedastic (p = 0.428).
This becomes more difficult with generalised models
Being able to interpret residual plots is important
Generalised Linear Models?
What are Generalised linear models?
The Generalized Linear Model (GLM) is an extension of the ordinary linear regression model that accommodates a broader range of response variable types and error distributions.
Key components:
Linear Predictor: Combines predictor variables linearly.
Link Function: Links the mean of the response variable to the linear predictor.
Error family: Drawn from the exponential family it determines the shape of the response variable’s distribution.
\(g\) the link function relating the linear predictor to the expected value (mean) of the response variable
\(\mu_i\) is the expected value of the response variable for observation \(i\)
\(f\) the inverse of the link function, transforming the linear predictor back to the response variable’s expected value
\(\epsilon_i\) is the error term for observation \(i\)
Changing model structures
There are a range of model structures available.
The choice of error family and link influence how well the model fits
Exponential Error Family
Link Function
Gaussian
Identity, sqrt, log, inverse
Binomial
Logit, probit, cloglog
Poisson
Log, identity, sqrt
Gamma
Inverse, identity, log, sqrt
Negative Binomial
Log, identity, sqrt
Questions?
A Linear Model Example
Janka timber hardness
This regression analysis uses data from the Australian forestry industry recording wood density and timber hardness from 36 different samples
Timber hardness is quantified on the “Janka” scale
The ‘amount of force required to embed a 0.444” steel ball into the wood to half of its diameter’.
Fitting the general linear model
janka_ls <-lm(hardness ~ dens, data = janka)summary(janka_ls)janka |>ggplot(aes( x = dens, y = hardness))+geom_point()+geom_smooth(method ="lm")
Call:
lm(formula = hardness ~ dens, data = janka)
Residuals:
Min 1Q Median 3Q Max
-338.40 -96.98 -15.71 92.71 625.06
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1160.500 108.580 -10.69 2.07e-12 ***
dens 57.507 2.279 25.24 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 183.1 on 34 degrees of freedom
Multiple R-squared: 0.9493, Adjusted R-squared: 0.9478
F-statistic: 637 on 1 and 34 DF, p-value: < 2.2e-16
Challenge
Fit the Linear Model lm(hardness ~ dens, data = janka) to the janka data.
Evaluate the model fit
10:00
Model diagnostics
plot(janka_ls, which=2)plot(janka_ls, which=3)
Violations of assumptions
library(lmtest)# Breusch-Pagan Test of Homoscedasticitylmtest::bptest(janka_ls)
studentized Breusch-Pagan test
data: janka_ls
BP = 4.4668, df = 1, p-value = 0.03456
# Shapiro-Wilk test for normality of residualsshapiro.test(residuals(janka_ls))
Shapiro-Wilk normality test
data: residuals(janka_ls)
W = 0.93641, p-value = 0.03933
Challenge
What can we try to improve the fit of our lm(hardness ~ dens, data = janka) ?
Let’s make suggestions then try them out
15:00
BoxCox
The R output for the MASS::boxcox() function plots a maximum likelihood curve (with a 95% confidence interval - drops down as dotted lines) for the best transformation for fitting the dependent variable to the model.
lambda value
transformation
0
log(Y)
0.5
sqrt(Y)
1
Y
2
Y^1
MASS::boxcox(janka_ls)
Square root transformation
janka_sqrt <-lm(sqrt(hardness) ~ dens, data = janka)plot(janka_sqrt, which=2)plot(janka_sqrt, which=3)
Model fit
ggplot(janka, aes(x = hardness, y = dens)) +geom_point() +# scatter plot of original data pointsgeom_smooth(method ="lm", formula = y ~sqrt(x)) +# regression linelabs(title ="Sqrt Linear Regression with ggplot2",x ="X", y ="Y") # axis labels and title
Log transformation
janka_log <-lm(log(hardness) ~ dens, data = janka)plot(janka_log, which=2)plot(janka_log, which=3)
Model fit
ggplot(janka, aes(x = hardness, y = dens)) +geom_point() +# scatter plot of original data pointsgeom_smooth(method ="lm", formula = (y ~log(x))) +# regression linelabs(title ="Log Linear Regression",x ="X", y ="Y") # axis labels and title
Polynomials
janka_poly <-lm(log(hardness) ~poly(dens, 2), data = janka)plot(janka_poly, which=2)plot(janka_poly, which=3)
Polynomial fit
summary(janka_poly)ggplot(janka, aes(x = hardness, y = dens)) +geom_point() +# scatter plot of original data pointsgeom_smooth(method ="lm", formula = (y ~poly(log(x), 2))) +# regression linelabs(title ="Quadratic Log Linear Regression",x ="X", y ="Y") # axis labels and title
Call:
lm(formula = log(hardness) ~ poly(dens, 2), data = janka)
Residuals:
Min 1Q Median 3Q Max
-0.22331 -0.05708 -0.01104 0.07500 0.18871
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7.1362 0.0168 424.896 < 2e-16 ***
poly(dens, 2)1 3.3862 0.1008 33.603 < 2e-16 ***
poly(dens, 2)2 -0.5395 0.1008 -5.354 6.49e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.1008 on 33 degrees of freedom
Multiple R-squared: 0.9723, Adjusted R-squared: 0.9706
F-statistic: 578.9 on 2 and 33 DF, p-value: < 2.2e-16
Weighted least squares
janka_wls <-lm(sqrt(hardness) ~ dens, weights =1/sqrt(hardness), data = janka)
Weights are specified as \(\frac{1}{\sqrt{\text{hardness}}}\), indicating that observations with higher hardness values will have lower weights, and vice versa.
This weighting scheme can assign more importance to certain observations in the model fitting process.
This approach is known as weighted least squares (WLS) regression, where observations are weighted differently based on certain criteria, such as the square root of hardness.
Call:
lm(formula = sqrt(hardness) ~ dens, data = janka, weights = 1/sqrt(hardness))
Weighted Residuals:
Min 1Q Median 3Q Max
-0.55072 -0.19009 -0.03393 0.19082 0.63845
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.01819 0.96647 2.088 0.0443 *
dens 0.76142 0.02201 34.601 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.294 on 34 degrees of freedom
Multiple R-squared: 0.9724, Adjusted R-squared: 0.9716
F-statistic: 1197 on 1 and 34 DF, p-value: < 2.2e-16
Troublesome transformations
Changes interpretability
Can introduce NEW violations of assumptions
Error and Mean change together
We don’t understand the true structure of our data
A Generalized linear model approach
Fit a glm(gaussian(link = "identity")) to the data
Check that the model is identical to the lm()
05:00
janka_glm <-glm(hardness ~ dens, data = janka, family =gaussian(link ="identity"))summary(janka_glm)
Call:
glm(formula = hardness ~ dens, family = gaussian(link = "identity"),
data = janka)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1160.500 108.580 -10.69 2.07e-12 ***
dens 57.507 2.279 25.24 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 33510.78)
Null deviance: 22485041 on 35 degrees of freedom
Residual deviance: 1139366 on 34 degrees of freedom
AIC: 481.21
Number of Fisher Scoring iterations: 2
A Generalized linear model approach
We previously identified that the square root of the dependent variable might have been a better fit with our linear model:
janka_glm <-glm(hardness ~ dens, data = janka, family =gaussian(link ="sqrt"))summary(janka_glm)
Call:
glm(formula = hardness ~ dens, family = gaussian(link = "sqrt"),
data = janka)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.60869 1.47660 1.767 0.0863 .
dens 0.75194 0.02722 27.622 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 25498.81)
Null deviance: 22485041 on 35 degrees of freedom
Residual deviance: 866960 on 34 degrees of freedom
AIC: 471.38
Number of Fisher Scoring iterations: 4
Plot model
ggplot(janka, aes(x = dens, y = hardness)) +geom_point() +# scatter plot of original data pointsgeom_smooth(method ="glm", method.args =list(gaussian(link ="sqrt"))) +# regression linelabs(title ="Linear Regression with ggplot2",x ="X", y ="Y") # axis labels and title
Residuals
janka_glm <-glm(hardness ~ dens, data = janka, family =gaussian(link ="sqrt"))plot(janka_glm, which=2)plot(janka_glm, which=3)
Fitting of the distribution ' gamma ' by maximum likelihood
Parameters :
estimate Std. Error
shape 3.351431971 0.4076228722
rate 0.002280106 0.0002166559
Loglikelihood: -287.9889 AIC: 579.9778 BIC: 583.1448
Correlation matrix:
shape rate
shape 1.0000000 0.7201126
rate 0.7201126 1.0000000
summary(fn)
Fitting of the distribution ' norm ' by maximum likelihood
Parameters :
estimate Std. Error
mean 1469.4722 131.71673
sd 790.3066 93.13779
Loglikelihood: -291.2889 AIC: 586.5779 BIC: 589.7449
Correlation matrix:
mean sd
mean 1 0
sd 0 1
Shape (α): Determines the shape or curve of the distribution. Higher values of α result in distributions that are more peaked and skewed to the left.
Rate (β): Influences the scale or spread of the distribution. Smaller values of β result in distributions that are more spread out, while larger values of β lead to distributions that are more concentrated around the mean.
Useful for modelling positively skewed continuous non-negative data
Gamma GLM
janka_gamma <-glm(hardness ~ dens, data = janka, family =Gamma(link ="sqrt"))summary(janka_gamma)
Call:
glm(formula = hardness ~ dens, family = Gamma(link = "sqrt"),
data = janka)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.86729 0.90357 2.067 0.0465 *
dens 0.76780 0.02243 34.224 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 0.009716208)
Null deviance: 11.26788 on 35 degrees of freedom
Residual deviance: 0.32876 on 34 degrees of freedom
AIC: 452.97
Number of Fisher Scoring iterations: 4
Plot model
ggplot(janka, aes(x = dens, y = hardness)) +geom_point() +# scatter plot of original data pointsgeom_smooth(method ="glm", method.args =list(Gamma(link ="sqrt"))) +# regression linelabs(title ="Linear Regression with ggplot2",x ="X", y ="Y") # axis labels and title
Exercise
Fit the Gamma distribution linear model to the data.
Evaluate the model fit against other glm fits.
10:00
Questions?
Likelihood
Log-Likelihood
Likelihood measures the probability of observing the given data under a specific statistical model, as a function of the model’s parameters.
It represents how well the model explains the observed data.
The log likelihood function in maximum likelihood estimations is usually computationally simpler as it is additive.
Although log-likelihood functions are mathematically easier than their multiplicative counterparts, they can be challenging to calculate by hand. They are usually calculated with software.
Maximum Likelihood Estimation
For a specified error family distribution and link function, the glm finds parameter values that have the highest likelihood score.
Which estimates have the are most probable for producing the observed data
This is an iterative process, without a single defined equation
For gaussian and identity link it will produce the same outcome as OLS.
The Akaike Information Criterion (AIC) balances the trade-off between model fit and complexity by penalizing models with more parameters.
It aims to find the model that best explains the data while avoiding overfitting.
AIC is calculated as follows: \(\text{AIC} = 2k - 2\ln(\hat{L})\), where: log-likelihood is the maximized value of the likelihood function for the fitted model. \(k\) is the number of parameters in the model.
In general a difference in AIC of <2 suggests equivalent goodness-of-fit
AIC
AIC can:
Compare non-nested models fitted to the same dataset and response variable
Compare across different GLM families and link functions
AIC cannot:
Compare models where the response variable is transformed
Compare models where there are differences in the dataset
Exercise
Revisit the janka models
Compare R-squared fits and AIC values
10:00
Likelihood ratio Test
The likelihood ratio test (LRT) “unifies” frequentist statistical tests. Brand-name tests like t-test, F-test, chi-squared-test, are specific cases (or even approximations) of the LRT.
The likelihood ratio test (LRT) begins with a comparison of the likelihood scores of the two models:
\(\text{LR} = -2(\ln L_1 - \ln L_2)\)
Where:
\(\ln L_1\) and \(\ln L_2\) are the log-likelihoods of the two models being compared.
Likelihood ratio test
This LRT statistic approximately follows a chi-square distribution.
To determine if the difference in likelihood scores among the two models is statistically significant, we next must consider the degrees of freedom.
In the LRT, degrees of freedom is equal to the number of additional parameters in the more complex model.
[1] 129.0546
Likelihood ratio test
The likelihood ratio test is typically calculated as twice the difference in deviance between the two models.
Scaling by -2 ensures the test statistic follows a chi-square distribution under certain conditions
We can calculate the likelihood ratio test statistic, compare it to the chi-square distribution with the appropriate degrees of freedom, and determine whether the improvement in fit between the two models is statistically significant.
Likelihood ratio test
Model 1: hardness ~ 1
Model 2: hardness ~ dens
#Df LogLik Df Chisq Pr(>Chisq)
1 2 -288.01
2 3 -223.48 1 129.05 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
F distribution
What if the variance is unknown e.g. Gaussian and Gamma distributions?
Then the null model estimates one parameter while the alternative model estimates two, so the difference in df is still 1.
In this case, we know this more familiarly as the ANOVA or F-test, which in this example is equivalent to the t-test on the intercept.
Call:
glm(formula = hardness ~ dens, family = Gamma(link = "sqrt"),
data = janka)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.86729 0.90357 2.067 0.0465 *
dens 0.76780 0.02243 34.224 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Gamma family taken to be 0.009716208)
Null deviance: 11.26788 on 35 degrees of freedom
Residual deviance: 0.32876 on 34 degrees of freedom
AIC: 452.97
Number of Fisher Scoring iterations: 4
Analysis of Deviance Table
Model 1: hardness ~ 1
Model 2: hardness ~ dens
Resid. Df Resid. Dev Df Deviance F Pr(>F)
1 35 11.2679
2 34 0.3288 1 10.939 1125.9 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Challenge?
Using the fruitfly.csv dataset
Build a model to investigate the effect of sleep, mating status and body size on longevity
Check model fits and report findings
30:00
Confidence intervals
# Assuming 'model' is your fitted GLM model# Extract coefficients and confidence intervalscoef_ci <-confint(fly_glm)# Extract estimatesestimates <-coef(fly_glm)# Combine estimates with confidence intervalsestimates_df <-data.frame(estimates = estimates,lower_ci = coef_ci[,1],upper_ci = coef_ci[,2])# Display dataframeestimates_df
Common response variable in ecological datasets is the binary variable: we observe a phenomenon X or its “absence”
Presence/Absence of a species
Presence/Absence of a disease
Success/Failure to observe behaviour
Survival/Death of organisms
Wish to determine if \(P/A∼ Variable\)
Called a logistic regression or logit model
Binary variables
Investigating mayfly abundance in relation to pollution gradient measured in Cumulative Criterion Units (CCU)
Mayflies serve as indicators of stream health due to sensitivity to metal pollution
Higher CCU values indicate increased metal pollution, potentially impacting mayfly presence or absence in the stream
Example
Fitting a linear regression vs a logit-link regression:
Bernoulli Distribution
The Bernoulli distribution models the probability of success or failure in a single trial of a binary experiment:
where success occurs with probability \(p\)
and failure with probability \({1-p}\)
Probability distribution
\(logit(p) = \log \frac{p}{1 - p}\)
\(y_i = Bernoulli(p) = \frac{p}{1 - p}\)
Mean of distribution Probability (p) of observing an outcome
Variance of observed responses As observed responses approach 0 or 1, the variance of the distribution decreases
The link function
If \(μ = xβ\) is only true for normally distributed data, then, if this is not the case, we must use a transformation on the expected values:
\[g(μ) = xβ\]
where \(g(μ)\) is the link function.
This allows us to relax the normality assumption.
The logit link
For binary data, the link function is called the logit:
\[ logit(p) = \log \frac{p}{1 - p} \]
The log of the odds \((p / (1 - p))\)
The logit link
\[ logit(p) = \log \frac{p}{1 - p} \]
The odds put our expected values on a 0 to +Inf scale.
The log transformation puts our expected values on a -Inf to +Inf scale.
Now, the expected values can be linearly related to the linear predictor.
Probability, odds, logit-odds
The odds put our expected values on a 0 to +Inf scale.
The log transformation puts our expected values on a -Inf to +Inf scale.
Now, the expected values can be linearly related to the linear predictor.
Example
mayfly_glm <-glm(Occupancy ~ CCU, family =binomial(link ="logit"), data = Mayflies)summary(mayfly_glm)
Call:
glm(formula = Occupancy ~ CCU, family = binomial(link = "logit"),
data = Mayflies)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.102 2.369 2.154 0.0313 *
CCU -3.051 1.211 -2.520 0.0117 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 34.795 on 29 degrees of freedom
Residual deviance: 12.649 on 28 degrees of freedom
AIC: 16.649
Number of Fisher Scoring iterations: 7
Each coefficient corresponds to a predictor variable, indicating the change in the log odds of the response variable associated with a one-unit change in that predictor, holding other predictors constant.
Interpret the model
Interpretation of coefficients involves exponentiating them to obtain odds ratios.
An odds ratio greater than 1 implies higher odds of the event occurring with an increase in the predictor.
An odds ratio less than 1 implies lower odds of the event occurring with an increase in the predictor.
Intepret the model
Probability
Odds
Log Odds
Verbiage
\(p=.95\)
\(\frac{95}{5} = \frac{19}{1} = 19\)
2.94
19 to 1 odds for
\(p=.75\)
\(\frac{75}{25} = \frac{3}{1} = 3\)
1.09
3 to 1 odds for
\(p=.50\)
\(\frac{50}{50} = \frac{1}{1} = 1\)
0
1 to 1 odds
\(p=.25\)
\(\frac{25}{75} = \frac{1}{3} = 0.33\)
-1.11
3 to 1 odds against
\(p=.05\)
\(\frac{95}{5} = \frac{1}{19} = 0.0526\)
-2.94
19 to 1 odds against
Logit odds
When we use a binomial model, we produce the ‘log-odds’, this produces a fully unconstrained linear regression as anything less than 1, can now occupy an infinite negative value -∞ to ∞.
We can interpret \(\beta_{1}\) and \(\beta_{2}\) as the increase in the log odds for every unit increase in \(x_{1}\) and \(x_{2}\).
We could alternatively interpret \(\beta_{1}\) and \(\beta_{2}\) using the notion that a one unit change in \(x_{1}\) as a percent change of \(e^{\beta_{1}}\) in the odds.
For the intercept:
The estimate as log odds is 5.102.
Therefore the odds are 164
Over 99% probability of observing mayfly when CCU = 0
Call:
glm(formula = Occupancy ~ CCU, family = binomial(link = "logit"),
data = Mayflies)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.102 2.369 2.154 0.0313 *
CCU -3.051 1.211 -2.520 0.0117 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 34.795 on 29 degrees of freedom
Residual deviance: 12.649 on 28 degrees of freedom
AIC: 16.649
Number of Fisher Scoring iterations: 7
For the predictor variable CCU:
The estimate is -3.051.
This means that for every one unit increase in CCU, the log odds of the response variable decreases by 3.051, holding all other predictors constant.
Call:
glm(formula = Occupancy ~ CCU, family = binomial(link = "logit"),
data = Mayflies)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 5.102 2.369 2.154 0.0313 *
CCU -3.051 1.211 -2.520 0.0117 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 34.795 on 29 degrees of freedom
Residual deviance: 12.649 on 28 degrees of freedom
AIC: 16.649
Number of Fisher Scoring iterations: 7
Now, let’s interpret the coefficient for CCU using the odds ratio interpretation:
The odds ratio associated with CCU is calculated as \(e^{\beta_{\text{CCU}}} = e^{-3.051}\).
This means that for every one unit increase in CCU, the odds of the response variable decrease by a multiple of 0.048
broom::tidy(mayfly_glm, exponentiate = T)
# A tibble: 2 × 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 164. 2.37 2.15 0.0313
2 CCU 0.0473 1.21 -2.52 0.0117
Challenge
Work with the odds and probability calculators
Check how comfortable you are with additive log-odds, multiplicative odds and changing probability
10:00
Model fit
For a simple Bernoulli/Binary GLM there are only a few checks that apply:
# Calculate and plot "deviance" residualsdresid <-resid(mayfly_glm, type ="deviance")hist(dresid)# Plot leverageplot(mayfly_glm, which =2)# Plot Cook's distanceplot(mayfly_glm, which =4)
Model fit
Residuals in binary glms will never be normally distributed - so a “residuals vs fitted plot” provides little insight.
DHARMa fit
DHARMa works by simulating residuals from the fitted model, comparing them with observed residuals, and providing diagnostic plots to assess model adequacy and identify potential issues like heteroscedasticity or influential data points.
Prediction
predict() function for GLMs computes predicted values based on the fitted model.
For existing data, it predicts values for each observation in the original dataset.
For new data, it predicts values for a specified new dataset, allowing for predictions on unseen data.
# Predictions for existing datapredict(mayfly_glm, type ="response")
# Predictions for new datanew_CCU<-data.frame(CCU =c(0,1,2,3,4,5))predict(mayfly_glm, newdata = new_CCU, type ="response")
Emmeans
The emmeans package in R computes estimated marginal means (EMMs) it can also provides post-hoc comparisons and contrasts for fitted models.
This specific function call below computes EMMs for the predictor variable CCU at specific values (0 to 5) and provides them on the response scale.
We can use emmeans package to compute estimated marginal means (EMMs) from a fitted GLM, convert them into a tibble format for visualization with ggplot2, overlay the means with uncertainty ribbons and add them on a scatter plot of the original data.
Exercise
Can you recreate the binomial glm for mayfly~ccu data?
10:00
Challenge
30:00
Using the malaria dataset fit a binomial glm to look at the predictors of malaria in the Seychelles warbler
Check the model fit
Make predictions
Produce a suitable figure
GLMs for Binomial data
Proportion data and GLM
Sometimes, count data aligns more closely with logistic regression methodologies than it initially appears.
We’re not looking at typical count data when measuring the number of occurrences with a known total sample size.
Imagine, for example, we’re studying the prevalence of a native underbrush species across various forest plots. We assess 15 different plots, each meticulously surveyed for the presence of this underbrush, counting the number of square meters where the underbrush is present versus absent within a standard area:
\[\frac{M^2~\text{with native underbrush (successes)}}{\text{Total surveyed area in}~M^2~(\text{trials})}\]
Bound between zero and one
Binomial distribution
\(logit(p) = \log \frac{p}{1 - p}\)
\(y_i = Binomial(n,p)\)
It is the collection of Bernoulli trials for the same event
It is represented by the number of Bernoulli trials \(n\)
It is also the probability of an event in each trial \(p\)
Binomial: Example
This small dataset is from an experiment looking at mortality in batches of flour exposed to different doses of pesticides.
The question of interest is:
How does pesticide dose affect mortality in the flour beetle Tribolium castaneum
Bernoulli deals with the outcome of the single trial of the event, whereas Binomial deals with the outcome of the multiple trials of the single event.
Bernoulli is used when the outcome of an event is required for only one time, whereas the Binomial is used when the outcome of an event is required multiple times.
Binomial glm
The GLM of the binomial counts will analyse the number of each batch of beetles killed, while taking into account the size of each group
Call:
glm(formula = cbind(dead, alive) ~ Dose, family = binomial(link = "logit"),
data = beetles)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -14.57806 1.29846 -11.23 <2e-16 ***
Dose 0.24554 0.02149 11.42 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 267.6624 on 7 degrees of freedom
Residual deviance: 8.4379 on 6 degrees of freedom
AIC: 38.613
Number of Fisher Scoring iterations: 4
Weights
An alternative method is to analyse the proportion of beetles killed and including the number of trials with the weight argument
Call:
glm(formula = Mortality_rate ~ Dose, family = binomial(link = "logit"),
data = beetles, weights = trials)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -14.57806 1.29846 -11.23 <2e-16 ***
Dose 0.24554 0.02149 11.42 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 267.6624 on 7 degrees of freedom
Residual deviance: 8.4379 on 6 degrees of freedom
AIC: 38.613
Number of Fisher Scoring iterations: 4
# Alternatively, using package-specific functions# Example with 'broom' packagelibrary(broom)tidy(beetle_glm, conf.int = T)
\(s\) corresponds to the number of times an event should occur.
\(p\) is the probability that the event will occur.
\((1 - p)\) is the probability that the event will not occur.
\(C\) term represents combinations, calculated as \(C(n, s) = \frac{n!}{s!(n - s)!}\), representing the number of ways to choose \(s\) successes out of \(n\) trials.
\(\theta\) term describes additional variance outside of the Binomial distribution
Fitting the quasibinomial model
Call:
glm(formula = cbind(dead, alive) ~ Dose, family = binomial(link = "logit"),
data = beetles)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -14.57806 1.29846 -11.23 <2e-16 ***
Dose 0.24554 0.02149 11.42 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 267.6624 on 7 degrees of freedom
Residual deviance: 8.4379 on 6 degrees of freedom
AIC: 38.613
Number of Fisher Scoring iterations: 4
Call:
glm(formula = cbind(dead, alive) ~ Dose, family = quasibinomial(link = "logit"),
data = beetles)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -14.57806 1.46611 -9.943 5.98e-05 ***
Dose 0.24554 0.02427 10.118 5.42e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for quasibinomial family taken to be 1.274895)
Null deviance: 267.6624 on 7 degrees of freedom
Residual deviance: 8.4379 on 6 degrees of freedom
AIC: NA
Number of Fisher Scoring iterations: 4
Fitting the quasibinomial model
Standard Error, Confidence intervals and p-values WILL change.
Estimates remain the same
Challenge
On the afternoon of January 28th, 1986 the space shuttle Challenger experienced a critical failure and broke apart in mid air, resulting in the deaths of all seven crewmembers.
An investigation discovered critical failures in all six of the O-rings in the solid rocket booster.
The question of interest is:
How does temperature affect the probability of o-ring failure?
What was the probability of all o-rings failing when the temperature was 36 degrees F
20:00
GLMs for count data
Poisson
Count or rate data are ubiquitous in the life sciences (e.g number of parasites per microlitre of blood, number of species counted in a particular area). These type of data are discrete and non-negative.
A useful distribution to model abundance data is the Poisson distribution:
a discrete distribution with a single parameter, λ (lambda), which defines both the mean and the variance of the distribution.
GLM
Recall the simple linear regression case (i.e a GLM with a Gaussian error structure and identity link). For the sake of clarity let’s consider a single explanatory variable where \(\mu\) is the mean for Y:
The mean function is unconstrained, i.e the value of \(\beta_0 + \beta_1X\) can range from \(-\infty\) to \(+\infty\).
Poisson
If we want to model count data we therefore want to constrain this mean to be positive only.
Mathematically we can do this by taking the logarithm of the mean (the log is the default link for the Poisson distribution).
We then assume our count data variance to be Poisson distributed (a discrete, non-negative distribution), to obtain our Poisson regression model (to be consistent with the statistics literature we will rename \(\mu\) to \(\lambda\)):
GLM Poisson
Poisson limitations
The Poisson distribution has the following characteristics:
Discrete variable (integer), defined on the range \(0, 1, \dots, \infty\).
A single rate parameter \(\lambda\), where \(\lambda > 0\).
Mean = \(\lambda\)
Variance = \(\lambda\)
Poisson: Case study
In a study by Kilner et al. (1999), the authors studied the begging rate of nestlings in relation to total mass of the brood of reed warbler chicks and cuckoo chicks. The question of interest is:
How does nestling mass affect begging rates between the different species?
This model is inadequate. It is predicting negative begging calls within the range of the observed data, which clearly does not make any sense.
cuckoo_lm <-lm(Beg ~ Mass + Species + Mass:Species, data = cuckoo)
Diagnostics
Let us display the model diagnostics plots for this linear model.
Curvature
Funnelling effect
Biological data and distributions
We can look at the kurtosis and skewness of our univariate variable - Begging:
We should therefore try a different model structure.
The response variable in this case is a classic count data: discrete and bounded below by zero (i.e we cannot have negative counts). We will therefore try a Poisson model using the canonical log link function for the mean:
λ varies with x (mass) which means residual variance will also vary with x, which means that we just relaxed the homogeneity of variance assumption!
Predicted values will now be integers instead of fractions
The model will never predict negative values (Poisson is strictly positive)
where \(M_i\) is nestling mass and \(S_i\) a dummy variable
\(S_i = \left\{\begin{array}{ll}1 & \text{if } i \text{ is warbler},\\0 & \text{otherwise}.\end{array}\right.\)
The term \(M_iS_i\) is an interaction term. Think of this as an additional explanatory variable in our model. Effectively it lets us have different slopes for different species (without an interaction term we assume that both species have the same slope for the relationship between begging rate and mass, and only the intercept differ).
Regression lines
The mean regression lines for the two species look like this:
cuckoo_glm1 <-glm(Beg ~ Mass + Species + Mass:Species, data=cuckoo, family=poisson(link="log"))summary(cuckoo_glm1)
Call:
glm(formula = Beg ~ Mass + Species + Mass:Species, family = poisson(link = "log"),
data = cuckoo)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.334475 0.227143 1.473 0.14088
Mass 0.094847 0.007261 13.062 < 2e-16 ***
SpeciesWarbler 0.674820 0.259217 2.603 0.00923 **
Mass:SpeciesWarbler -0.021673 0.008389 -2.584 0.00978 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 970.08 on 50 degrees of freedom
Residual deviance: 436.05 on 47 degrees of freedom
AIC: 615.83
Number of Fisher Scoring iterations: 6
Note there appears to be a negative interaction effect for Species:Mass, indicating that Begging calls do not increase with mass as much as you would expect for Warblers as compared to Cuckoos.
Exercise
Fit a Poisson model to the cuckoo data
Look at the residual plots - how have they changed?
10:00
Model summary
summary(cuckoo_glm1)
Call:
glm(formula = Beg ~ Mass + Species + Mass:Species, family = poisson(link = "log"),
data = cuckoo)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.334475 0.227143 1.473 0.14088
Mass 0.094847 0.007261 13.062 < 2e-16 ***
SpeciesWarbler 0.674820 0.259217 2.603 0.00923 **
Mass:SpeciesWarbler -0.021673 0.008389 -2.584 0.00978 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 970.08 on 50 degrees of freedom
Residual deviance: 436.05 on 47 degrees of freedom
AIC: 615.83
Number of Fisher Scoring iterations: 6
VIF, or Variance Inflation Factor, measures the extent of multicollinearity among independent variables in regression analysis.
High multicollinearity, indicated by elevated VIF values, can distort regression coefficients and inflate standard errors, undermining the reliability of the model’s predictions
car::vif(cuckoo_glm1)
Mass Species Mass:Species
4.050091 12.967267 14.496767
check_model(cuckoo_glm1, check ="vif")
Mean-centering
True VIF is highly problematic
Mean centering can remove apparent VIF as means become zero
In Poisson GLMs, using an offset allows us to model rates rather than counts, by including the logarithm of the exposure variable as an offset, ensuring that the model accounts for the differing exposure levels across observations.
Call:
glm(formula = Call_Total ~ mass.c + Species + mass.c:Species,
family = poisson(link = "log"), data = cuckoo, offset = log(Mins))
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.249441 0.030831 72.961 <2e-16 ***
mass.c 0.129261 0.002502 51.671 <2e-16 ***
SpeciesWarbler 0.300202 0.035330 8.497 <2e-16 ***
mass.c:SpeciesWarbler -0.057893 0.002829 -20.462 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 9597.9 on 50 degrees of freedom
Residual deviance: 3071.8 on 47 degrees of freedom
AIC: 3334.6
Number of Fisher Scoring iterations: 6
Poisson limitations
So for the Poisson regression case we assume that the mean and variance are the same. Hence, as the mean increases, the variance increases also (heteroscedascity). This may or may not be a sensible assumption so watch out! Just because a Poisson distribution usually fits well for count data, doesn’t mean that a Gaussian distribution can’t always work.
Recall the link function between the predictors and the mean and the rules of logarithms (if \(\log{\lambda} = k\)(value of predictors), then \(\lambda = e^k\)):
Overdispersion
When the residual deviance is higher than the residual degrees of freedom, we say that the model is overdispersed. This situation is mathematically represented as:
\[
\phi = \frac{\text{Residual deviance}}{\text{Residual degrees of freedom}}
\]
Overdispersion occurs when the variance in the data is even higher than the mean, indicating that the Poisson distribution might not be the best choice. This can be due to many reasons, such as the presence of many zeros, many very high values, or missing covariates in the model.
Overdispersion and what to do about it
Overdispersion statistic values > 1
Causes of over-dispersion
What to do about it
Model mis-specification (missing covariates or interactions)
Add more covariates or interaction terms
Too many zeros (“zero inflation”)
Use a zero-inflated model
Non-independence of observations
Use a generalised mixed model with random effects to take non-independence into account
Variance is larger than the mean
Use a quasipoisson GLM if overdispersion = 2-15. Use a negative binomial GLM if > 15-20
Overdispersion and what to do about it
It may be true that you are initially unsure whether overdispersion is coming from zero-inflation, larger variance than the mean or both.
You can fit models with the same dependent and independent variables and different error structures and use AIC to help determine the best model for your data
Overdispersion in the cuckoo model
summary(cuckoo_glm2)
Call:
glm(formula = Beg ~ mass.c + Species + mass.c:Species, family = poisson(link = "log"),
data = cuckoo)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.263284 0.091946 24.615 < 2e-16 ***
mass.c 0.094847 0.007261 13.062 < 2e-16 ***
SpeciesWarbler 0.234068 0.106710 2.194 0.02827 *
mass.c:SpeciesWarbler -0.021673 0.008389 -2.584 0.00978 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 970.08 on 50 degrees of freedom
Residual deviance: 436.05 on 47 degrees of freedom
AIC: 615.83
Number of Fisher Scoring iterations: 6
Quasi-Poisson GLM
The systematic part and the link function remain the same
\[
\log{\lambda} = \beta_0 + \beta_1 x_i
\]
phi\(\phi\) is a dispersion parameter, it will not affect estimates but will affect significance. Standard errors are multiplied by \(\sqrt\phi\)
\[
Y_i = Poisson(\phi \lambda_i)
\]
Where:
\(Yi\) is the response variable.
\(\phi\) is the dispersion parameter, which adjusts the variance of the distribution
Confidence intervals and p-values WILL change.
Exercise
Fit a quasipoisson model to the cuckoo data
cuckoo_glm2 <-glm(Beg ~ mass.c + Species + mass.c:Species, data=cuckoo, family=quasipoisson(link="log"))
Look at the residual plots - how have they changed?
Look at the SE and p-values - how have they changed
Calculate new 95% confidence intervals
F distributions
The presence of overdispersion suggested the use of the F-test for nested models. We will test if the interaction term can be dropped from the model.
Though here for a term with only one degree of freedom \(t^2\) = the F distribution with denominator degree of freedom
drop1(cuckoo_glm2, test ="F")
Single term deletions
Model:
Beg ~ mass.c + Species + mass.c:Species
Df Deviance F value Pr(>F)
<none> 436.05
mass.c:Species 1 442.81 0.7293 0.3974
Zero-inflation
“Zero-inflation” refers to the problem of “too many zeros”.
The dependent variable contains more zeros than would be expected under the standard statistical distribution for zero-bounded count data (Poisson or negative binomial).
Zero-inflation is a common feature of count datasets. Large numbers of zeros can arise because
There was nothing there to count (true zeros)
Things were there but missed (false zeros)
Zero-inflation
Model
Acronym/Abbr
Source of Overdispersion
Zero-inflated Poisson
ZIP
Zeros
Negative Binomial
NegBin
Positive integers
Zero-inflated Negative Binomial
ZINB
Zeros + Positive integers
ZIP in R
In a ZIP model we define two halves of the model; the first half models the count process (i.e. a Poisson GLM explaining non-zero integers) and the second half models the binomial (i.e. explains the zeros).
cuckoo_zip <-zeroinfl(Beg ~ mass.c + Species + mass.c:Species| mass.c + Species + mass.c:Species,dist ="poisson",link ="logit",data = cuckoo)summary(cuckoo_zip)
Call:
zeroinfl(formula = Beg ~ mass.c + Species + mass.c:Species | mass.c +
Species + mass.c:Species, data = cuckoo, dist = "poisson", link = "logit")
Pearson residuals:
Min 1Q Median 3Q Max
-3.2611 -0.9164 -0.1979 0.8878 3.7917
Count model coefficients (poisson with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.395981 0.104443 32.515 < 2e-16 ***
mass.c 0.027481 0.008805 3.121 0.001802 **
SpeciesWarbler -0.800911 0.117996 -6.788 1.14e-11 ***
mass.c:SpeciesWarbler 0.038046 0.009819 3.875 0.000107 ***
Zero-inflation model coefficients (binomial with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.8893 0.7431 1.197 0.2314
mass.c -0.2154 0.1078 -1.998 0.0457 *
SpeciesWarbler -8.6156 4.1849 -2.059 0.0395 *
mass.c:SpeciesWarbler -0.3357 0.3483 -0.964 0.3352
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Number of iterations in BFGS optimization: 18
Log-likelihood: -155.6 on 8 Df
Model validation
Model validation is not straightforward for zero-inflated models
sresid <-residuals(cuckoo_zip, type ="pearson")pred <-predict(cuckoo_zip)
par(mfrow =c(2, 2), mar =c(4, 4, 2, 1))hist(sresid)plot(sresid ~ pred)qqnorm(sresid)qqline(sresid, col ="red")
Negative Binomial
Negative Binomial
Negative binomial GLMs are favored when overdispersion is high.
It has two parameters, \(μ\) and \(\theta\). \(\theta\) controls for the dispersion parameter (smaller \(\theta\) indicates higher dispersion). It corresponds to a combination of two distributions (Poisson and gamma).
It assumes that the \(Y_i\) are Poisson distributed with the mean \(μ_i\) assumed to follow a gamma distribution:
Negative binomial is not in the glm() function. We need the MASS package:
Suitable for strongly over-dispersed zero-bounded integer data.
Negbin: Example
The example is an experiment measuring the effect of the parasitic tapeworm Schistocephalus solidus infection on the susceptibility of infection from a second parasite, the trematode Diplostomum pseudospathaceum
The question of interest is:
How does Treatment/Infection status affect trematode counts?
Exercise
10:00
Determine whether the poisson model is suitable for this dataset?
